# 德劳内三角网实现算法，支持各种几何
# 前置： pip3 install shapely,pyshp
import shapefile, geojson,math,json
from shapely.ops import triangulate,unary_union,split
from shapely import wkt
from shapely.geometry import MultiPoint,multipolygon,Polygon,MultiPolygon,shape
import shapely.geometry as geometry
import numpy as np

def Delaunay_simple():
    basePath = r'C://Users//srhee//Desktop//poly//poly2//'
    path = basePath + 'poly2.shp'
    sf = shapefile.Reader(path)
    shapes = sf.shapes()
    w = shapefile.Writer(basePath + 'polygonSimple')
    w.field('name', 'C')

    for i in range(len(shapes)):
        shp = sf.shape(i)
        poly = Polygon(shp.points)
        print(poly)

        triangles = triangulate(poly, tolerance=0.00001, edges=False)
        for j in triangles:
            g1 = wkt.loads(j.wkt)
            g2 = geojson.Feature(geometry=g1, properties={})
            # print(g2.geometry.coordinates)
            w.poly(g2.geometry.coordinates)
            w.record('polygon')
    w.close()

def Delaunay_holes():
    basePath = r'C://Users//srhee//Desktop//poly//polygon1//'
    path = basePath + 'polygon1.shp'
    sf = shapefile.Reader(path)
    shapes = sf.shapes()
    w = shapefile.Writer(basePath + 'polygonHoles')
    w.field('name', 'C')

    for i in range(len(shapes)):
        shp = sf.shape(i)

        points = Polygon(shp.points)
        # print(points)
        # g12 = geojson.Feature(geometry=wkt.loads(points.wkt), properties={})
        g12 = geojson.Feature(geometry=shp, properties={})
        shellCoor = g12.geometry.coordinates[:1]
        holesCoor = g12.geometry.coordinates[1:]
        # print(shellCoor)
        # print(holesCoor)

        shellTuple = tuple(tuple([y for y in x]) for x in shellCoor[0])
        holesTupleArr = []
        for m in holesCoor:
            holesTuple = tuple(tuple([y for y in x]) for x in m)
            holesTupleArr.append(holesTuple)

        # print(shellTuple)

        # triangulate 不支持带洞polygon进行三角剖分
        # 变相进行判断，总减去洞，这里生成的三角面含所有顶点
        trianglesAll = triangulate(Polygon(shellTuple, holesTupleArr), tolerance=0.00001, edges=False)
        # trianglesAll = triangulate(points, tolerance=0.00001, edges=False)
        trianglesHoles = []
        for m in holesCoor:
            holesTuple = tuple(tuple([y for y in x]) for x in m)
            trianglesHoles = trianglesHoles + triangulate(Polygon(holesTuple), tolerance=0.00001, edges=False)

        # 坐标一致则重复，标识并记录总的集合下标
        indexArr = []
        for j in trianglesAll:
            geo_j = geojson.Feature(geometry=wkt.loads(j.wkt), properties={})

            for k in trianglesHoles:
                geo_k = geojson.Feature(geometry=wkt.loads(k.wkt), properties={})
                if geo_j.geometry.coordinates == geo_k.geometry.coordinates:
                    indexArr.append(trianglesAll.index(j))

        # 这里注意要把下标集合进行从大到小进行排序，再del删除
        for y in sorted(indexArr,reverse=True):
            del trianglesAll[y]

        for j in trianglesAll:
            g1 = wkt.loads(j.wkt)
            print('图形wkt：',g1,'面积：',j.area)
            g2 = geojson.Feature(geometry=g1, properties={})
            w.poly(g2.geometry.coordinates)
            # w.shapeType = shapefile.POLYGON
            w.record('polygon')

    w.close()

def Delaunay_simple_v2():
    basePath = r'C://Users//srhee//Desktop//poly//poly2//'
    path = basePath + 'poly2.shp'
    sf = shapefile.Reader(path)
    shapes = sf.shapes()
    w = shapefile.Writer(basePath + 'polygonSimple2')
    w.field('name', 'C')

    for i in range(len(shapes)):
        shp = sf.shape(i)
        poly = Polygon(shp.points)
        print(poly)

        triangles = triangulate(poly, tolerance=0.00001, edges=False)

        temp = []
        for j in triangles:
            # 加一个判断 中心点
            if poly.intersects(j.centroid):
                temp.append(j)
                # g1 = wkt.loads(j.wkt)
                # g2 = geojson.Feature(geometry=g1, properties={})
                # # print(g2.geometry.coordinates)
                # w.poly(g2.geometry.coordinates)
                # w.record('polygon')
        print( unary_union(temp) )
        print(  )
        test = poly.difference(unary_union(temp))
        # print( split( poly ,unary_union(temp)).wkt )
        g2 = geojson.Feature(geometry = wkt.loads(test.wkt), properties={})
        # print(g2.geometry.coordinates)
        w.poly(g2.geometry.coordinates)
        w.record('polygon')

    w.close()

def Delaunay_holes_v2():
    basePath = r'C://Users//srhee//Desktop//poly//polygon1//'
    path = basePath + 'polygon1.shp'
    sf = shapefile.Reader(path)
    shapes = sf.shapes()
    w = shapefile.Writer(basePath + 'polygonHoles')
    w.field('name', 'C')

    for i in range(len(shapes)):
        shp = sf.shape(i)

        points = Polygon(shp.points)
        # print(points)
        # g12 = geojson.Feature(geometry=wkt.loads(points.wkt), properties={})
        g12 = geojson.Feature(geometry=shp, properties={})
        shellCoor = g12.geometry.coordinates[:1]
        holesCoor = g12.geometry.coordinates[1:]
        # print(shellCoor)
        # print(holesCoor)

        shellTuple = tuple(tuple([y for y in x]) for x in shellCoor[0])
        holesTupleArr = []
        for m in holesCoor:
            holesTuple = tuple(tuple([y for y in x]) for x in m)
            holesTupleArr.append(holesTuple)

        # print(shellTuple)

        # triangulate 不支持带洞polygon进行三角剖分
        # 变相进行判断，总减去洞，这里生成的三角面含所有顶点
        trianglesAll = triangulate(Polygon(shellTuple, holesTupleArr), tolerance=0.00001, edges=False)
        # trianglesAll = triangulate(points, tolerance=0.00001, edges=False)
        trianglesHoles = []
        for m in holesCoor:
            holesTuple = tuple(tuple([y for y in x]) for x in m)
            trianglesHoles = trianglesHoles + triangulate(Polygon(holesTuple), tolerance=0.00001, edges=False)

        # 坐标一致则重复，标识并记录总的集合下标
        indexArr = []
        for j in trianglesAll:
            geo_j = geojson.Feature(geometry=wkt.loads(j.wkt), properties={})

            for k in trianglesHoles:
                geo_k = geojson.Feature(geometry=wkt.loads(k.wkt), properties={})
                if geo_j.geometry.coordinates == geo_k.geometry.coordinates:
                    indexArr.append(trianglesAll.index(j))

        # 这里注意要把下标集合进行从大到小进行排序，再del删除
        for y in sorted(indexArr,reverse=True):
            del trianglesAll[y]

        for j in trianglesAll:
            g1 = wkt.loads(j.wkt)
            print('图形wkt：',g1,'面积：',j.area)
            g2 = geojson.Feature(geometry=g1, properties={})
            w.poly(g2.geometry.coordinates)
            # w.shapeType = shapefile.POLYGON
            w.record('polygon')

    w.close()

if __name__ == '__main__':
    # 生成简单图形三角面
    # Delaunay_simple()
    # 生成带洞图形三角面
    # Delaunay_holes()

    # 2.0 升级版，支持凹多边形
    # 生成简单图形三角面
    Delaunay_simple_v2()
    # 生成带洞图形三角面
    # Delaunay_holes()
